### Replication data for: Ethnic Minorities, Interstate War, and Popular Support for Fiscal Capacity Development
### Author: Andre Walter & Patrick Emmenegger
### Contact: andre.walter[at]ipz.uzh.ch

### Install rdrobust package version 2.0.1 to replicate the results (package "devtools" required)

devtools::install_version("rdrobust", version = "2.0.1", repos = "http://cran.us.r-project.org")

### Packages

library(rdrobust)
library(tidyverse)
library(sf)
library(modelsummary)
library(cowplot)
library(xtable)
library(ggrepel)

rm(list=ls())

### Data

load("findat.RData")

load("mapdata.RData")

load("mapadd1.RData")

load("mapadd2.RData")

##### Figure 1: RDD Plots

ggplot(findat, aes(x = dist, y = elev, color = factor(north))) +
  geom_point(size = 1, alpha = 0.5) + 
  geom_smooth(data = filter(findat, dist < 0), method = "lm",formula=y ~ poly(x, 3, raw=TRUE)) +
  geom_smooth(data = filter(findat, dist >= 0), method = "lm",formula=y ~ poly(x, 3, raw=TRUE)) +
  geom_vline(xintercept = 0) +
  labs(x = "Distance from Natural Border", y = "Elevation") + theme_bw()+theme(legend.position = "none")


ggplot(findat, aes(x = dist, y = yes_sh, color = factor(north))) +
  geom_point(size = 1, alpha = 0.5) + 
  geom_smooth(data = filter(findat, dist < 0), method = "lm",formula=y ~ poly(x, 3, raw=TRUE)) +
  geom_smooth(data = filter(findat, dist >= 0), method = "lm",formula=y ~ poly(x, 3, raw=TRUE)) +
  geom_vline(xintercept = 0) +
  labs(x = "Distance from Natural Border", y = "Yes Share") + theme_bw()+theme(legend.position = "none")


ggplot(findat, aes(x = dist, y = yes_sh2, color = factor(north))) +
  geom_point(size = 1, alpha = 0.5) + 
  geom_smooth(data = filter(findat, dist < 0), method = "lm",formula=y ~ poly(x, 3, raw=TRUE)) +
  geom_smooth(data = filter(findat, dist >= 0), method = "lm",formula=y ~ poly(x, 3, raw=TRUE)) +
  geom_vline(xintercept = 0) +
  labs(x = "Distance from Natural Border", y = "Yes Share") + theme_bw()+theme(legend.position = "none")


ggplot(findat, aes(x = dist, y = fra_sh, color = factor(north))) +
  geom_point(size = 1, alpha = 0.5) + 
  geom_smooth(data = filter(findat, dist < 0), method = "lm",formula=y ~ poly(x, 3, raw=TRUE)) +
  geom_smooth(data = filter(findat, dist >= 0), method = "lm",formula=y ~ poly(x, 3, raw=TRUE)) +
  geom_vline(xintercept = 0) +
  labs(x = "Distance from Natural Border", y = "Share French Speaker") + theme_bw()+
  theme(legend.position = "none")


##### Figure 2: Maps of Bern with yes share direct taxation 1915/1918

### Vote 1915

ggplot() +
  geom_sf(data = BE, aes(fill = yes_sh2)) + theme_map()+
  labs(fill = "Yes Share Direct Taxes")+
  ggsn::scalebar(BE, dist_unit = "km", dist = 25, st.size=3, height=0.01, transform = T)+
  #geom_sf(data = border, color = "#CC6600", size = 1.5)+
  geom_sf(data = border, color = "black", size = 1)+
  scale_fill_distiller(palette = "Oranges", type="seq", direction = 1)+
  geom_text_repel(data = border, aes(x = X,y = Y,label = Name), 
                  fontface = "bold", nudge_x = c(1), nudge_y = c(.21))+
  geom_sf(data = jurabord,color = "black", size = 1)+
  geom_text_repel(data = jurabord, aes(x = X,y = Y,label = Name), 
                  fontface = "bold", nudge_x = c(.15), nudge_y = c(1))#+ theme(plot.margin = margin(0,-30, 0, -30, "cm"))


# Vote 1918

ggplot() +
  geom_sf(data = BE, aes(fill = yes_sh)) + theme_map()+
  labs(fill = "Yes Share Direct Taxes")+
  ggsn::scalebar(BE, dist_unit = "km", dist = 25, st.size=3, height=0.01, transform = T)+
  #geom_sf(data = border, color = "#CC6600", size = 1.5)+
  geom_sf(data = border, color = "black", size = 1)+
  scale_fill_distiller(palette = "Oranges", type="seq", direction = 1)+
  geom_text_repel(data = border, aes(x = X,y = Y,label = Name), 
                  fontface = "bold", nudge_x = c(1), nudge_y = c(.21))+
  geom_sf(data = jurabord,color = "black", size = 1)+
  geom_text_repel(data = jurabord, aes(x = X,y = Y,label = Name), 
                  fontface = "bold", nudge_x = c(.15), nudge_y = c(1))#+ theme(plot.margin = margin(0,-20, 0, -30, "cm"))



##### Table 1: Covariate imbalance statistics

varnam <- c()
est <- c()
pval <- c()
rdd <- c()

test <- rdrobust(findat$yes_sh, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$yes_sh2, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$elev, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$fra_sh, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$agri, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$cowagri, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$cath_sh, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$popdens, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$poor_sh, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$Grund_mean, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$Naehr_mean, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$Skelet_mea, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$Wasser_mea, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$Wassersp_m, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$Vernaess_m, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$rough_mean, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$Wheat_mean, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$Potatoes_m, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))
test <- rdrobust(findat$distbe, findat$dist)
est <- c(est, test[[1]][[1]])
pval <- c(pval,test[[6]][[3]])
rdd <- c(rdd,"summary")
varnam <- c(varnam,as.character(test$call$y[[3]]))

test <- data.frame(code = varnam,est = est,pval = round(pval,3), adpval = round(p.adjust(pval, method = "fdr"),3))

test <- data.frame(rdd = rdd, code = varnam,est = est,pval = round(pval,3), adpval = round(p.adjust(pval, method = "fdr"),3))

myres <- test %>% filter(rdd == "summary") %>% select(est, adpval) %>% 
   rename(Estimate = est, p.value = adpval) 
  
bastat <- findat %>% mutate(treat = ifelse(dist>0,1,0)) 

Nobs <- function(x) length(x)-sum(is.na(x))

options("modelsummary_format_numeric_latex" = "plain")

datasummary(Heading("Support Direct Taxation 1918")*yes_sh+Heading("Support Direct Taxation 1915")*yes_sh2+Heading("Elevation")*elev+
              Heading("Share French Speaker")*fra_sh + Heading("Share Agricultural Employment")*agri+
              Heading("ln Cows per Agricultural Worker")*cowagri +
              Heading("Share Catholics")*cath_sh + 
              Heading("Population Density")*popdens+
              Heading("Share Poor")*poor_sh+Heading("Soil Thoroughness")*Grund_mean +
              Heading("Nutrient Storage Capacity")*Naehr_mean+ Heading("Stone Content")*Skelet_mea+
              Heading("Water Permeability")*Wasser_mea+
              Heading("Water Storage Capacity")*Wassersp_m+Heading("Water Logging")*Vernaess_m+
              Heading("Terrain Roughness")*rough_mean+Heading("Climate Wheat")*Wheat_mean+
              Heading("Climate Potatoes")*Potatoes_m+Heading("Distance Capital")*distbe  ~ 
              factor(treat)*((Mean+Min+Max+SD)*Arguments(na.rm=TRUE,fmt = "%.3f")), bastat, 
            title = "Covariate Imbalance Statistics", add_columns = myres)

##### Figure 3: Main estimation

### Vote 1915

n <- 20
m <- 6

resh <- matrix(NA, ncol = m, nrow = n)

refu <- matrix(NA, ncol = m, nrow = n)

for (i in c(4:n)){
  
  sharp <- rdrobust(findat$yes_sh2, findat$dist, h = i)
  
  resh[i,] <- c(0,i, sharp$Estimate[[1]], sharp$ci[[1]], sharp$ci[[4]],sum(sharp$N_h))
  
  fuzzy <- rdrobust(findat$yes_sh2, findat$dist, h = i, fuzzy = findat$fra_sh)
  
  refu[i,] <- c(1,i, fuzzy$Estimate[[1]], fuzzy$ci[[1]], fuzzy$ci[[4]],sum(fuzzy$N_h))
  
}

rddres <- bind_rows(data.frame(resh),data.frame(refu))

rddres <- rddres %>% drop_na() %>% rename(fuzzy = X1, bw = X2, est = X3, low = X4, upp = X5, obs = X6) %>%
  mutate(bwop = ifelse(bw == round(rdrobust(findat$yes_sh, findat$dist)$bws[[1]],0) & fuzzy == 0, 1,
                       ifelse(bw == round(rdrobust(findat$yes_sh, findat$dist, fuzzy = findat$fra_sh)$bws[[1]],0) & fuzzy == 1, 1,0)),
         bwop = factor(bwop), fuzzy = ifelse(fuzzy == 1, "Fuzzy", "Sharp"))



estfig <- ggplot(data = rddres, aes(y = est, x = bw, group = fuzzy, linetype = fuzzy, color = as.factor(bwop)))+
  geom_point(aes(),position=position_dodge(0.5))+theme_bw()+
  geom_errorbar(aes(ymin = low, ymax = upp),position=position_dodge(0.5))+scale_color_manual(values = c("black","grey"),guide = FALSE)+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())+labs(title = "", x = "", group = "RDD Type",
                                             y = "Treatment Effect")+ geom_hline(yintercept = 0, size = 0.2)+
  scale_linetype_manual(name = "RDD Type", values = c(1,2))

estobs <- ggplot(data = rddres, aes(y = obs, x = bw, group = fuzzy, linetype = fuzzy))+geom_line(aes())+theme_bw()+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())+labs(title = "", x = "Bandwidth",
                                             y = "Number of Observations")+ geom_hline(yintercept = 0, size = 0.2)+
  scale_linetype_manual(name = "RDD Type", values = c(1,2))

plot_grid(estfig,estobs, nrow = 2, align = "v", rel_heights = c(3, 1))


### Vote 1918

n <- 20
m <- 6

resh <- matrix(NA, ncol = m, nrow = n)

refu <- matrix(NA, ncol = m, nrow = n)

for (i in c(4:n)){

sharp <- rdrobust(findat$yes_sh, findat$dist, h = i)

resh[i,] <- c(0,i, sharp$Estimate[[1]], sharp$ci[[1]], sharp$ci[[4]],sum(sharp$N_h))

fuzzy <- rdrobust(findat$yes_sh, findat$dist, h = i, fuzzy = findat$fra_sh)

refu[i,] <- c(1,i, fuzzy$Estimate[[1]], fuzzy$ci[[1]], fuzzy$ci[[4]],sum(fuzzy$N_h))
  
}

rddres <- bind_rows(data.frame(resh),data.frame(refu))

rddres <- rddres %>% drop_na() %>% rename(fuzzy = X1, bw = X2, est = X3, low = X4, upp = X5, obs = X6) %>%
  mutate(bwop = ifelse(bw == round(rdrobust(findat$yes_sh, findat$dist)$bws[[1]],0) & fuzzy == 0, 1,
                ifelse(bw == round(rdrobust(findat$yes_sh, findat$dist, fuzzy = findat$fra_sh)$bws[[1]],0) & fuzzy == 1, 1,0)),
         bwop = factor(bwop), fuzzy = ifelse(fuzzy == 1, "Fuzzy", "Sharp"))


estfig <- ggplot(data = rddres, aes(y = est, x = bw, group = fuzzy, linetype = fuzzy, color = as.factor(bwop)))+
  geom_point(aes(),position=position_dodge(0.5))+theme_bw()+
  geom_errorbar(aes(ymin = low, ymax = upp),position=position_dodge(0.5))+scale_color_manual(values = c("black","grey"),guide = FALSE)+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())+labs(title = "", x = "", group = "RDD Type",
                                             y = "Treatment Effect")+ geom_hline(yintercept = 0, size = 0.2)+
  scale_linetype_manual(name = "RDD Type", values = c(1,2))

estobs <- ggplot(data = rddres, aes(y = obs, x = bw, group = fuzzy, linetype = fuzzy))+geom_line(aes())+theme_bw()+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())+labs(title = "", x = "Bandwidth",
                                             y = "Number of Observations")+ geom_hline(yintercept = 0, size = 0.2)+
  scale_linetype_manual(name = "RDD Type", values = c(1,2))

plot_grid(estfig,estobs, nrow = 2, align = "v", rel_heights = c(3, 1))

##### Table A1: Main estimation

rddtab <- rddres %>% select(-bwop) %>% group_split(fuzzy) %>% bind_cols() %>% select(-`fuzzy...1`,-`fuzzy...7`)

colnames(rddtab) <- c("Bandwidth","Point Estimate","Lower CI","Upper CI","Num. Obs.","Bandwidth","Point Estimate","Lower CI","Upper CI","Num. Obs.")

tabx <- xtable(rddtab, caption = "Main RDD Estimates", label = "rddtab", 
               align = c("c","c","c","c","c","c","c","c","c","c","c"),digits=c(0,0,3,3,3,0,0,3,3,3,0))
print(tabx, booktabs = TRUE, hline.after = c(-1,0,nrow(tabx)), include.rownames = F,
      floating.environment = "sidewaystable")


##### Figure 4: Placebo Jura

n <- 20
m <- 5

resh <- matrix(NA, ncol = m, nrow = n)


for (i in c(3:n)){
  
  sharp <- rdrobust(findat$yes_sh, findat$juradist, h = i)
  
  resh[i,] <- c(i, sharp$Estimate[[1]], sharp$ci[[1]], sharp$ci[[4]],sum(sharp$N_h))
  
}

rddres <- bind_rows(data.frame(resh))

rddres <- rddres %>% drop_na() %>% rename(bw = X1, est = X2, low = X3, upp = X4, obs = X5) %>%
  mutate(bwop = ifelse(bw == round(rdrobust(findat$yes_sh, findat$juradist)$bws[[1]],0) , 1,0),
         bwop = factor(bwop))

estfig <- ggplot(data = rddres, aes(y = est, x = bw, color = as.factor(bwop)))+
  geom_point()+theme_bw()+
  geom_errorbar(aes(ymin = low, ymax = upp),)+scale_color_manual(values = c("black","grey"),guide = FALSE)+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())+labs(title = "", x = "", y = "Treatment Effect")+ 
  geom_hline(yintercept = 0, size = 0.2)

estobs <- ggplot(data = rddres, aes(y = obs, x = bw))+geom_line(aes())+theme_bw()+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())+labs(title = "", x = "Bandwidth",
                                             y = "Number of Observations")+ geom_hline(yintercept = 0, size = 0.2)

plot_grid(estfig,estobs, nrow = 2, align = "v", rel_heights = c(3, 1))


##### Figure A2:: RDD with controls (1918, share of the poor excluded due to many missings)

covar <- cbind(findat$agri,findat$cath_sh,findat$popdens,findat$Grund_mean,findat$Naehr_mean, 
               findat$rough_mean,findat$Wheat_mean,findat$Potatoes_m,findat$Skelet_mea,
               findat$Wasser_mea,findat$Vernaess_m,findat$Wassersp_m,findat$cowagri)

n <- 20
m <- 6

resh <- matrix(NA, ncol = m, nrow = n)

refu <- matrix(NA, ncol = m, nrow = n)

for (i in c(4:n)){
  
  sharp <- rdrobust(findat$yes_sh, findat$dist, h = i, covs = covar)
  
  resh[i,] <- c(0,i, sharp$Estimate[[1]], sharp$ci[[1]], sharp$ci[[4]],sum(sharp$N_h))
  
  fuzzy <- rdrobust(findat$yes_sh, findat$dist, h = i, fuzzy = findat$fra_sh, covs = covar)
  
  refu[i,] <- c(1,i, fuzzy$Estimate[[1]], fuzzy$ci[[1]], fuzzy$ci[[4]],sum(fuzzy$N_h))
  
}

rddres <- bind_rows(data.frame(resh),data.frame(refu))

rddres <- rddres %>% drop_na() %>% rename(fuzzy = X1, bw = X2, est = X3, low = X4, upp = X5, obs = X6) %>%
  mutate(bwop = ifelse(bw == round(rdrobust(findat$yes_sh, findat$dist)$bws[[1]],0) & fuzzy == 0, 1,
                       ifelse(bw == round(rdrobust(findat$yes_sh, findat$dist, fuzzy = findat$fra_sh)$bws[[1]],0) & fuzzy == 1, 1,0)),
         bwop = factor(bwop), fuzzy = ifelse(fuzzy == 1, "Fuzzy", "Sharp"))



estfig <- ggplot(data = rddres, aes(y = est, x = bw, group = fuzzy, linetype = fuzzy, color = as.factor(bwop)))+
  geom_point(aes(),position=position_dodge(0.5))+theme_bw()+
  geom_errorbar(aes(ymin = low, ymax = upp),position=position_dodge(0.5))+scale_color_manual(values = c("black","red"),guide = FALSE)+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())+labs(title = "", x = "", group = "RDD Type",
                                             y = "Treatment Effect")+ geom_hline(yintercept = 0, size = 0.2)+
  scale_linetype_manual(name = "RDD Type", values = c(1,2))

estobs <- ggplot(data = rddres, aes(y = obs, x = bw, group = fuzzy, linetype = fuzzy))+geom_line(aes())+theme_bw()+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())+labs(title = "", x = "Bandwidth",
                                             y = "Number of Observations")+ geom_hline(yintercept = 0, size = 0.2)+
  scale_linetype_manual(name = "RDD Type", values = c(1,2))

plot_grid(estfig,estobs, nrow = 2, align = "v", rel_heights = c(3, 1))

##### Figure A3: Alternative bandwidth estimators

m <- 6

bwest <- c("cerrd","certwo","cersum","cercomb1","cercomb2")

n <- length(bwest)

resh <- matrix(NA, ncol = m, nrow = n)

refu <- matrix(NA, ncol = m, nrow = n)

resh2 <- matrix(NA, ncol = m, nrow = n)

refu2 <- matrix(NA, ncol = m, nrow = n)

for (i in 1:length(bwest)){
  
  sharp <- rdrobust(findat$yes_sh, findat$dist, bwselect = bwest[i])
  
  resh[i,] <- c("Sharp Linear",bwest[i], sharp$Estimate[[1]], sharp$ci[[1]], sharp$ci[[4]],sum(sharp$N_h))
  
  fuzzy <- rdrobust(findat$yes_sh, findat$dist, fuzzy = findat$fra_sh, bwselect = bwest[i])
  
  refu[i,] <- c("Fuzzy Linear",bwest[i], fuzzy$Estimate[[1]], fuzzy$ci[[1]], fuzzy$ci[[4]],sum(fuzzy$N_h))
  
  sharp2 <- rdrobust(findat$yes_sh, findat$dist, bwselect = bwest[i], p = 2)
  
  resh2[i,] <- c("Sharp Quadratic",bwest[i], sharp2$Estimate[[1]], sharp2$ci[[1]], sharp2$ci[[4]],sum(sharp2$N_h))
  
  fuzzy2 <- rdrobust(findat$yes_sh, findat$dist, fuzzy = findat$fra_sh, bwselect = bwest[i], p = 2)
  
  refu2[i,] <- c("Fuzzy Quadratic",bwest[i], fuzzy2$Estimate[[1]], fuzzy2$ci[[1]], fuzzy2$ci[[4]],sum(fuzzy2$N_h))
  
}

results <- bind_rows(data.frame(resh),data.frame(refu),data.frame(resh2),data.frame(refu2))


ggplot(results, aes(x = X2, y = as.numeric(X3), linetype = factor(X1), group = factor(X1))) + 
  geom_point(position=position_dodge(0.5)) + geom_errorbar(aes(ymin = as.numeric(X4), ymax = as.numeric(X5)), width = 0.2,position=position_dodge(0.5)) +
  geom_hline(yintercept = 0) + labs(x = "Bandwidth Estimator", y = "Treatment Effect")+
  scale_linetype_manual(name = "Estimator", values = c(1:4))+theme_bw()+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())

##### Figure A4: RDD arbitrary placebos

varnam <- c()
est <- c()
lci <- c()
uci <- c()

above <- findat %>% filter(dist>0)

test <- rdrobust(above$yes_sh, above$dist, c = quantile(above$dist,.25))
est <- c(est, test$Estimate[[1]])
lci <- c(lci,test$ci[[1]])
uci <- c(uci,test$ci[[4]])
varnam <- c(varnam,"North First Quartile")

test <- rdrobust(above$yes_sh, above$dist, c = median(above$dist))
est <- c(est, test$Estimate[[1]])
lci <- c(lci,test$ci[[1]])
uci <- c(uci,test$ci[[4]])
varnam <- c(varnam,"North Median")

test <- rdrobust(above$yes_sh, above$dist, c = quantile(above$dist,.75))
est <- c(est, test$Estimate[[1]])
lci <- c(lci,test$ci[[1]])
uci <- c(uci,test$ci[[4]])
varnam <- c(varnam,"North Third Quartile")


below <- findat %>% filter(dist<0)

test <- rdrobust(below$yes_sh, below$dist, c = quantile(below$dist,.25))
est <- c(est, test$Estimate[[1]])
lci <- c(lci,test$ci[[1]])
uci <- c(uci,test$ci[[4]])
varnam <- c(varnam,"South First Quartile")

test <- rdrobust(below$yes_sh, below$dist, c = median(below$dist))
est <- c(est, test$Estimate[[1]])
lci <- c(lci,test$ci[[1]])
uci <- c(uci,test$ci[[4]])
varnam <- c(varnam,"South Median")

test <- rdrobust(below$yes_sh, below$dist, c = quantile(below$dist,.75))
est <- c(est, test$Estimate[[1]])
lci <- c(lci,test$ci[[1]])
uci <- c(uci,test$ci[[4]])
varnam <- c(varnam,"South Third Quartile")

test <- data.frame(code = varnam, est = est, lci = lci, uci = uci)

ggplot(test, aes(x = code, y = est)) + 
  geom_point() + geom_errorbar(aes(ymin = lci, ymax = uci), width = 0.2) +
  geom_hline(yintercept = 0) + labs(x = "Placebo Cutoff", y = "Treatment Effect")+
  theme_bw()+
  theme(axis.line.x = element_line(colour = "black"), axis.line.y = element_line(colour = "black"),
        panel.border = element_blank())

